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ABSTRACT 


An  experimental  and  analytical  investigation  was  conducted  to  develop  a  calibrated 
torsional  model  of  a  three-cylinder,  two-stroke  diesel  engine.  A  Detroit  Diesel  3-53  engine 
was  instrumented  for  time  resolved  measurement  of  cylinder  firing  pressures  and  high 
resolution  near  instantaneous  shaft  speed  using  a  720  and  a  3,600  count  per  revolution  optical 
encoder.  Data  were  taken  for  three  speeds  and  three  torques  for  a  total  of  nine  conditions.  A 
six  degree-of-freedom  torsional  vibration  model  of  the  crankshaft,  connecting  rods,  and 
pistons  was  developed.  The  nonlinear  inertias,  due  to  the  reciprocating  pistons,  were  included 
along  with  linear  stiffness  and  damping.  The  equations  of  motion  were  numerically  integrated 
over  a  cycle  to  obtain  predicted  response.  The  predicted  response  was  compared  to  the 
measured  response  at  the  free  end  of  the  crankshaft. 
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I.  INTRODUCTION 

A.         MOTIVATION 

Diesel  engines  are  widely  used  in  both  military  and  civilian  applications.  In  the 
Navy,  they  are  utilized  for  main  propulsion  and  electric  power  generation.  The  reliable 
detection  of  faults  both  non-obtrusively  and  inexpensively  could  shift  the  focus  of 
maintenance  from  regularly  scheduled  to  condition  based.  This  maintenance  strategy  has 
the  potential  to  save  costly,  unneeded  maintenance  and  can  increase  availability  of  diesel 
engines. 

Individual  cylinder  pressures  of  internal  combustion  engines  are  excellent 
predictors  of  many  types  of  engine  faults.  Cylinder  pressures  are  directly  related  to  the 
torque  output  capability  of  an  engine.  Therefore,  the  contribution  of  torque  from  an 
individual  cylinder  is  directly  related  to  the  maximum  pressure  that  cylinder  develops. 
Cylinders  experiencing  loss  of  pressure  can  be  caused  by  faulty  fuel  injectors,  worn 
cylinder  lining,  faulty  piston  rings,  leaking  valves,  and  blown  or  leaking  cylinder  head 
gaskets. 

There  is  a  problem  with  measuring  cylinder  firing  pressures  in  the  field.  While  high 
quality  pressure  transducers  can  accurately  measure  cylinder  pressures,  they  have  a  limited 
lifetime  in  the  combustion  chamber  environment  and  are  expensive.  Many  alternative 
indirect  methods  of  cylinder  pressure  monitoring  such  as  the  use  of  strain  bolts  are  being 
explored. 


Another  indirect  method  is  to  measure  the  variations  in  shaft  speed.  With  this 
information  and  equations  of  motion  describing  the  dynamics  of  the  engine,  the  firing 
pressures  are  solved  for  numerically.  In  order  to  accurately  predict  cylinder  firing 
pressures  with  shaft  speed  variation,  an  accurate  torsional  engine  model  representative  of 
the  physical  system  is  required. 

B.         STATE  OF  THE  ART 

There  have  been  numerous  research  efforts  focused  on  extracting  cylinder  pressure 
information  from  crankshaft  speed  fluctuations.  It  is  well  known  that  these  speed 
fluctuations  are  a  direct  result  of  the  cyclic  nature  of  the  crankshaft  torque  waveform. 
This  torque  waveform  is  a  direct  result  from  gas  pressure  forces,  due  to  combustion,  and 
inertia  forces,  due  to  internal  reciprocating  masses.  The  challenge  has  been  taking  these 
measured  speed  fluctuations  and  passing  them  through  a  robust  engine  model  to  predict 
cylinder  pressures.  This  literature  review  focused  in  areas  of  non-obtrusive  cylinder 
pressure  and  torque  prediction  with  emphasis  on  torsional  vibration  modeling. 

Wilson  [Ref.  1]  developed  a  torsional  vibration  model  of  a  crankshaft  with 
reciprocating  components.  With  his  methods,  the  individual  inertia,  stiffness,  and  damping 
between  cylinder  crankthrows  and  the  rest  of  the  shaft  are  readily  determined.  His 
methods  are  currently  used  by  industry  to  develop  torsional  vibration  models  of 
crankshafts.  Kabele  [Ref.  2]  developed  a  torsional  vibration  model  of  a  diesel  crankshaft 
which  included  piston  and  ring  friction,  connecting  rod  plane  motion,  hysteretic  damping 
in  the  crankshaft,  and  cylinder  pressure  as  a  forcing  function.    The  model  was  derived 


from  Newtonian  Mechanics  and  the  governing  differential  equations  were  non-linear 
making  a  numerical  solution  necessary.  Once  his  model  was  developed,  he  compared 
measured  and  model  predicted  vibration  amplitudes  for  various  harmonic  orders  and 
concluded  that  correlation  between  the  two  was  excellent.  The  measured  vibration 
amplitudes  were  from  a  V-8  diesel  engine. 

Citron  et  al,  [Ref  3]  did  a  computer  simulation  and  produced  cylinder  pressure 
torque  waveforms  by  using  an  elastic  model  of  an  engine,  drivetrain  system,  and  measured 
speed  fluctuations.  These  speed  fluctuation  data  were  passed  back  through  the  model  to 
determine  the  fluctuation  waveform  of  both  the  total  engine  torque  being  developed  and 
the  cylinder  pressure  waveform  which  gave  rise  to  it. 

Rizzoni  [Ref.  4]  proposed  passing  crankshaft  speed  fluctuations  through  an 
equivalent  electronic  circuit  representing  engine  dynamics  to  predict  cylinder  misfires.  He 
conducted  experimental  work  on  production  spark  ignited  vehicles  and  confirmed  that  it  is 
possible  to  apply  this  method  with  minimal  hardware  and  computational  overhead. 

Sobel  et  al,  [Ref.  5]  described  a  technique  for  measuring  instantaneous  crankshaft 
torque  utilizing  a  non-contact  ferromagnetic  material  to  measure  torsional  stress  in  the 
crankshaft.  The  instantaneous  torque  of  an  internal  combustion,  four-cylinder  engine  was 
measured  at  the  flywheel  and  compared  with  pressure  signals  during  bench  tests 
Excellent  correlation  between  this  non-contact  measured  torque  and  values  of  the  mean 
effective  pressure  for  each  cylinder  were  obtained. 

Brown  and  Neill  [Ref  6]  presented  a  non-contact  method  for  determining 
simultaneously  the  pressure  in  each  cylinder  by  employing  a  pattern  recognition  technique 


to  compare  crankshaft  speed  fluctuations  to  reference  patterns  in  a  knowledge  base.  The 
non-contact  method  utilized  was  an  interval  timer  and  a  magnetic  sensor  which  timed  the 
flywheel  gear  teeth  as  they  passed.  The  experiment  was  carried  out  on  a  Detroit  Diesel 
6V-92TA  engine.  The  experimental  results  show  that  the  method  evaluates  the  cylinder 
pressures  with  an  RMS  error  of  less  than  six  percent.  But  this  requires  extensive  data 
bases  for  each  condition. 

Bell  [Ref.  7]  concluded,  utilizing  the  same  Detroit  Diesel  3-53  engine  used  in  this 
research,  that  the  variation  of  shaft  speed  held  information  which  could  be  used  to  predict 
cylinder  firing  pressures.  He  also  recommended  the  development  of  a  torsional  vibration 
model  of  the  engine  and  flywheel  system  along  with  a  more  rigid  optical  encoder  mounting 
for  more  accurate  time  resolved  measurement. 

C.         METHODOLOGY  AND  OBJECTIVES 

The  engine  used  for  this  research  is  a  Detroit  Diesel  Series  53  engine,  model  5033- 
500 IN.  It  is  a  three-cylinder,  two-stroke  engine.  During  one  revolution,  the  crankshaft 
speed  fluctuates  as  a  result  of  cylinder  pressure  variations.  These  crankshaft  speed 
fluctuations  can  be  measured  using  an  optical  encoder.  A  torsional  vibration  model  with 
nonlinear  inertias,  due  to  the  reciprocating  pistons,  and  linear  stiffness  and  damping  will  be 
developed  to  represent  the  dynamics  of  the  engine.  The  torsional  vibration  model  will  be 
calibrated  by  ensuring,  as  closely  as  possible,  predicted  response  matches  the  measured 
response.  Data  from  these  measured  speed  fluctuations  can  then  later  be  passed  through 
the  calibrated  torsional  model  to  determine  cylinder  firing  pressures 


The  objectives  of  this  thesis  are  to: 

1 .  Instrument  a  diesel  engine  with  a  rigid  optical  encoder  mount  such  that  reliable  time 
resolved  measurements  of  crankshaft  speed  relative  to  the  engine  block  can  be  obtained. 

2.  Develop  a  torsional  vibration  model  of  the  crankshaft  including  components  connected 
to  it  (i.e.  connecting  rods,  pistons,  flywheels,  and  external  loads). 

3.  Operate  the  diesel  engine  at  several  speeds  and  applied  torques  to  establish  baseline 
data. 

4.  Compare  the  model  predicted  response  to  the  experimental  measured  response. 

D.         ORGANIZATION 

Chapter  II  describes  the  development  of  the  crankshaft  torsional  vibration  model 
Equations  of  motion  are  derived  for  a  six  degree-of-freedom  torsional  model  and 
numerical  solution  procedure  is  presented. 

Chapter  III  describes  the  engine,  dynomometer,  instrumentation,  and  data 
acquisition  procedures. 

Chapter  IV  presents  results  of  data  acquisition,  the  measured  response,  and  the 
model  predicted  response.  A  comparison  of  measured  to  predicted  response  is  made 
based  on  the  results. 

Chapter  V  contains  summary,  conclusions,  and  recommendations  for  future  work 


II.  TORSIONAL  VIBRATION  MODEL 

A.         PHYSICAL  SYSTEM 

The  engine  used  for  this  research  is  representative  of  diesel  engines.  It  is  a  three- 
cylinder,  two-stroke  engine.  Engine  operating  conditions  and  design  specifications  are 
discussed  in  more  detail  in  Section  III.  A 

Figure  2.1  illustrates  the  engine  layout  with  the  crankshaft,  flywheel,  pistons  and 
connecting  rods.  The  front  of  the  engine  is  to  the  left  The  cylinders  are  numbered  from 
front  (1)  to  back  (3).  Connected  to  the  power  take-off  shaft  on  the  right  is  the  load, 
which  in  this  case  is  supplied  by  a  water  brake  dynamometer.  The  firing  order  is  1-3-2  and 
the  crankshaft  rotates  clockwise  as  seen  from  the  front 
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Figure  2. 1    Detroit  Diesel  3-53  with  crankshaft  and  major  reciprocating 
components.  The  front  of  the  engine  is  the  left  hand  side. 
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The  crankshaft  crankthrows  are  the  locations  where  the  connecting  rods  are  attached 
to  the  crankshaft.  Crankthrows  number  one,  two,  and  three  are  displaced  120  degrees. 
After  the  number  one  piston  reaches  top  dead  center  (TDC),  its  highest  point  of  travel, 
120  degrees  of  crankshaft  rotation  later  the  number  three  piston  arrives  at  its  TDC. 
Number  two  piston's  TDC  is  240  degrees  after  number  one.  The  power  into  the 
crankshaft  is  delivered  at  each  crankthrow. 

The  forces  due  to  gas  pressure  in  each  cylinder  can  be  taken  as  a  force  '-F'  on  the 
piston  and  an  equal  and  opposite  force  'F'  on  the  cylinder  head.  The  force  '-F'  is 
transmitted  through  the  connecting  rod  to  the  crankthrow.  The  resultant  force  on  the 
crankshaft  multiplied  by  the  eccentricity  'R'  of  the  crankthrow  creates  the  gas  torque.  If 
there  were  no  friction,  the  gas  torque  would  equal  the  shaft  applied  torque.  Figure  2.2 
shows  the  relationship  between  gas  pressure  force  and  torque. 
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where:  F  =  cylinder  pressure  *  area  of  the  piston  (lbf) 
R  =  crankthrow  eccentricity  (in) 
0  =  angular  position  of  crankshaft  relative  to  TDC  (rad) 
L  =  length  of  connecting  rod  (in) 


Figure  2.2  Relationship  between  gas  pressure  forces  and  gas  torque  into  the  crankshaft. 


Gas  torque  into  the  crankshaft  as  a  result  of  cylinder  gas  pressures  is  given  by  the 
following  relationship: 


Ticyi=PicyiRsin(^)Ap  where:  P#cyi  =  Measured  cylinder  pressure  (psia) 

T2cyi=P2cyiRsin(  0  -  47t/3)Ap  R  =  Crankshaft  crankthrow  eccentricity  (in) 

T3Cyi=P3cyiRsin(#-27r./3)Ap  6  =  Angular  position  of  #1  crankthrow 

after  TDC  (rad) 
(2.1)  Ap  =  Cross-sectional  area  of  the  pistons  (in2) 


The  torque  applied  to  the  crankshaft  is  somewhat  less  than  the  gas  torques  due 
mainly  to  friction  on  the  cylinder  wall  (rings).  The  major  components  considered  in  the 
development  of  a  torsional  vibration  model  are  the  crankshaft,  pistons,  connecting  rods, 
and  flywheel.  This  physical  system  to  be  modeled  is  shown  in  Figure  2.3. 
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Figure  2.3  Crankshaft  with  a  single  reciprocating  component. 
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1.  Model  Development 

The  reciprocating  components  of  the  engine  are  very  complex.  However,  to 
facilitate  calculation  it  is  useful  to  represent  the  complex  system  by  a  reduced  order 
system.  For  this  research,  six  rigid  masses  connected  by  sections  of  shafting  with  torsional 
stiflhess  and  damping  shall  be  used.  This  reduced  system  should  retain  the  dominant 
dynamic  characteristics  of  the  original  physical  system  that  are  of  interest.  The  six  degree- 
of-freedom  torsional  vibration  model  that  is  used  to  represent  the  physical  system  is  shown 
in  Figure  2.4. 
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Figure  2.4  Torsional  vibration  model  of  crankshaft  and  reciprocating  components 


In  Figure  2.4,  0i  represents  the  angular  position  of  the  crankshaft  where  the  optical 
encoder  attaches,  02  through  04  represents  the  angular  position  of  crankthrows  number 
one  through  three,  05  represents  the  angular  position  of  the  flywheel,  and  06  represents  the 
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angular  position  of  the  load  (water  brake  dynamometer).  Each  'J'  represents  the 
equivalent  mass  polar  moment  of  inertia  at  each  degree  of  freedom.  Equivalent  torsional 
stiffness  between  each  degree  of  freedom  is  represented  by  Ki  through  K5.  The  material 
damping  of  the  crankshaft  between  each  degree  of  freedom  is  represented  by  C12,  C23,  C34, 
C45,  and  C56.  The  linear  viscous  damping  effects  of  each  cylinder  due  to  translation  of  the 
pistons  in  the  cylinders  and  oil  at  the  journal  bearing  surfaces  are  represented  by  C2,  C3, 
and  C4.  All  inertias  (J's)  and  stiffnesses  (K's)  are  readily  determined  by  the  geometry  and 
material  properties  of  the  system.  Procedures  for  determining  the  equivalent  inertia  and 
stiffiiess  for  the  model  are  outlined  in  Appendix  A.  The  damping  for  the  system,  especially 
viscous  damping  effects,  are  not  so  readily  determined  and  will  be  discussed  later  when 
comparing  measured  results  to  model  predicted  results.  Table  2. 1  defines  the  values  of 
inertia  (J's)  and  stiffness  (K's)  for  the  model  [Ref  8]. 

Table  2. 1  Inertia  and  Stiffness  of  the  model. 


INERTIA 

STIFFNESS 

J,  =  0.029  lbrin-sec2 

K,  -  3.470  x  106  lbrin/rad 

h  =  0.273  +  O.O69[l-cos(20)J  lbrin-sec2 

K2  =  9.551  x  106  lbpin/rad 

J3  =  0.171  +  0.069J  l-cos2(9-27t/3)]  lbrin-sec2 

K3  =  9.551  x  106  lbrin/rad 

J4  =  0.273  +  O.O69[1-cos2(0-4ti/3)]  lbrin-sec2 

JL,  =  13.49  x  106  lbrin/rad 

J5  =  6.194  lbrin-sec2 

K5=  1.304  xlO6  lbrin/rad 

J6  =  0.280  lbrin-sec2 

2.  Data  Required  For  Model 

Physical  dimensions  were  obtained  by  measuring  the  geometry  of  an  equivalent 
crankshaft  out  of  an  engine  of  the  same  model  as  shown  in  Figure  2.3.  The  material  used 
to  manufacture  the  crankshaft  is  Steel  SAE  1548  Modified  (Maganese  1.00/1.30  and 
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Carbon  0.43/0.52).       Both  the  physical  dimensions  and  material  properties  were  checked 
against  the  manufacturer's  design  drawings  [Ref.  9] 

Other  data  required  for  the  model  are  measured  gas  torque  as  a  function  of  angular 
position  and  load  torque 

B.         EQUATIONS  OF  MOTION 

Referring  to  Figure  2.4,  the  equations  governing  the  motion  of  the  system  are: 

(1)  Ji0i  +  Ci2(0i-02)  +  £i(0i-02)  =  O 

(2)  {J  2  +  J2rgc  )02  +  Cn(02  -  01 )  +  ATl(6>2  -  0l)  +  Cll(0l  -  03)  +  Kl(0l  -  03)  +  ClOl  =  Tlcyl  (t) 

(3)  (J 3  +  Jirec  )03  +  C23(03  -  02)  +  /C2(03  -  02)  +  C34(03  -  04)  +  £3(03  -  04)  +  C303  =  Tlcyl  (t) 

(4)  (J 4  +  J4rec  )04  +  C34(04  -  03)  +  Ki(04  ~  03)  +  C45(04  -  05)  +  AT<04  -  05)  +  C404  =  Tscyl  (t) 

(5)  J505  +  C45(05  -  04)  +  £4(05  -  04)  +  C56(05  -  06)  +  £5(05  -  06>  =  0 

(6)  J 606  +  C56(06  -  05)  +  Ki(06  -  05)  =  -Tload  (2.2) 

The  equations  of  motion  in  matrix  format  are: 
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The  reciprocating  component  of  inertia  at  each  crankthrow,  J^ec,  J3rec,  and  J4rec, 
represent  the  inertia  contribution  of  the  piston  and  connecting  rod.  The  influence  of  these 
reciprocating  components  is  a  maximum  when  the  crankshaft  is  at  the  position 
corresponding  to  mid-stroke  (90  degrees  after  top  dead  center),  and  disappears  when  the 
crank  is  on  top  or  bottom  dead  center  (0  or  180  degrees  after  top  dead  center).  The 
relationship  of  the  non-linear  reciprocating  inertia  to  crankshaft  position  is: 


J2rec=  WR2/2g[  1  -cos(2  6 )]       where:  W  =  Weight  of  piston  and  connecting  rod  ( 1 0. 5  lbf) 


J3rec=  WR2/2g[l-cos2(0-27t/3)] 
J4rec=WR2/2g[  1  -cos2(  6  -4tc/3  ] 

(2.4) 


R  =  Crankshaft  crankthrow  eccentricity  (2.25  in) 
g  =  Acceleration  of  gravity  (386  in/sec2) 
0  =  Angular  position  of  #1  crankthrow  after 
TDC  (rad) 
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The  right  hand  side  of  the  equations  of  motion  are  the  applied  torques  [Ticyi(t), 
T2cyi(t),  and  T3cyi(t)]  into  the  system  and  load  torque,  Tioad  (which  is  considered  constant), 
out  of  the  system.  Figure  2.5  shows  an  example  of  individual  measured  cylinder  gas 
torque  contribution  referenced  to  top  dead  center  of  number  one  cylinder.  Since  the  total 
gas  torque  into  the  system  has  to  overcome  damping  and  friction  before  it  exits  as  useful 
torque  (load  torque),  the  gas  torques  can  be  multiplied  by  some  factor  to  account  for  the 
losses.  It  is  also  possible  to  subtract  a  constant  parasitic  torque  from  each  cylinder  gas 
torque.  Integrating  each  cylinder  torque  curve  with  respect  to  angular  position  yields  and 
average  torque.  Dividing  the  load  torque  by  the  sum  of  the  three  average  cylinder  torques 
will  yield  the  mechanical  efficiency  (rjm)  to  be  used.  Specifically: 


Tin=^-JT(t)d0      and     nm=  =    T-f    -       (2.5) 
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Figure  2.5  Individual  cylinder  gas  torque  contribution  referenced  to  Nr.  1 
cylinder  top  dead  center  (TDC). 


C.         NUMERICAL  SOLUTION  OF  EQUATIONS  OF  MOTION 

In  order  to  solve  the  six  second  order  equations  of  motion  it  is  first  necessary  to 
reduce  them  to  a  system  of  twelve  first  order  differential  equations.  These  are: 
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o7 

-ox 

o* 

=  K 

o9 

--oz 

^10 

=  0* 

On 

-o5 

0n 

=  'o6 

Jx07  -  C]2  0%  +  CX207  -  Kx02  +  KX0X  -  0 

V \  +  J2nc  W%  -  c23e9  +  (c23  +  cX2  +  c2  )6S  -  cl2e7  -  K2e3  +  (kx  +  k2  )$2  -kxox  =  TXcyl(t) 

(/j  +  ^3rec)09  ~  Q4^10  +  (Q4  +  C23  +  Q  )^9  ~  Qb^S  _  ^3^4  +  (^2  +  ^3^3  "  K2^2   =  ^cy/CO 

W ;  +^Ao  -QA  +(c„  +c34  +  c4)^10  -qa  -^5  +(^3  +^4)^4  -^3  =  r3^(o 

J50n  -  C560X2  +  (C56  +  C45 )0XX  -  C450,o  -K506  +  (K4  +  K5 )0S  -  K404  =  0 

J<P\1  +  ^56^12  -  Q6^ll  +  K5&6  -  K5^5   =  ~Tload 

(2.6) 

This  initial  value  problem  was  solved  using  a  fourth  order  Runge-Kutta  method 
ordinary  differential  equation  solver  (ODE45  in  MATLAB).  Twelve  initial  conditions  are 
required,  namely,  angular  position  0{t  =  0)  and  angular  velocity  0{t  =  0) .  The  system 
repeats  after  one  revolution  of  the  crankshaft  (2k  rad).  A  numerical  solution  for  the 
predicted  angular  position  of  the  first  degree  of  freedom  is  calculated  and  presented  in 
Section  IV.  B. 
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III.  EXPERIMENTAL  SETUP 

A.         THE  THREE  CYLINDER  TWO  STROKE  DIESEL  ENGINE 
1.  Operating  Conditions  and  Design  Specifications 

The  test  engine  used  is  a  Detroit  Diesel  Series  53  two  stroke  three  cylinder  engine. 
The  engine  design  operating  characteristics  are  listed  in  Table  3.1. 

Table  3.1  Engine  Characteristics  [Ref.  10] 


Model 

5033-5001N 

Engine  Type 

In  line-2  Cycle-Naturally  Aspirated 

Number  of  Cylinders 

3 

Bore  and  Stroke 

3.875x4.50  inches 

Exhaust  Valves  per  cylinder 

4 

Engine  Displacement 

159  cubic  inches 

Compression  Ratio 

21.0:1 

Maximum  Power  Output 

92bhp 

Full  Load  Speed 

2,800  RPM 

Peak  Torque 

198  ft-lbs 

Peak  Torque  Speed 

1,500  RPM 

Brake  Mean  Effective  Pressure 

o->    lb 

83  — 
in 
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The    manufacturer    provides    the    following    estimate    of   cylinder    combustion 
pressures  based  upon  an  empirically  obtained  relationship: 


r       4,650,000  IHP  ... 

C  „  =  — = where:  CD  =  Combustion  pressure  (psia), 

p      B2SNRPM  P  F  F 


(3.1) 


N  =  No.  of  cylinders 

B  =  Bore  (in),  S  =  Stroke  (in), 

RHP 

IHP  =  indicated  horsepower  ( =  t]mech  =  0.70 ) 

IHP 

RPM  =  revolutions  per  minute  of  crankshaft 
This  formula  gives  values  within  ±  5%  for  speeds  ( 1 200  RPM  to  maximum  rated 
RPM)  and   power  (continuous  to   maximum   intermittent   rated   power).      [Ref.    11] 
Appendix  B  [Ref.  12]  provides  rated  power  as  a  function  of  speed  curves  and  engine 
specification  data. 

2.  Engine  Event  Timing 

The  major  events  occurring  in  the  cylinders  with  respect  to  crankshaft  angle  are 
given  in  Figure  3.1,  specifically,  intake  port,  exhaust  valve,  and  fuel  injector  movements 
with  respect  to  top  dead  center  of  the  number  one  cylinder.  The  firing  order  of  the  engine 
(1-3-2)  can  clearly  be  seen.  The  number  three  cylinder  and  number  one  cylinder  events 
occur  120  and  240  degrees  respectively  after  the  corresponding  number  one  cylinder 
events 


#3  Cyl  Intake  Pons 

#3  Cyl  Exhaust 
#3  Cyl  Injector 

#2  Cyl  Intake  Ports 

#2  Cyl  Exhaust 

#2  Cyl  Iniector 

#1  Cyl  Intake  Ports 

#1  Cyl  Exhaust 

#1  Cyl  Iniector 

180        150        120        90         60         30         C  30         SO         90         120        ISO        180 


Figure  3.1   Engine  event  timing  map  referenced  to  TDC  Nr.  1  cylinder 
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B  INSTRUMENTATION 

1.  Optical  Encoder  and  Mount 

The  optical  encoder  mount  holds  the  casing  of  the  optical  encoder  while  the  shaft 
of  the  encoder  is  connected  to  the  crankshaft.  The  mount  is  in  turn  attached  to  the  engine 
block.  Any  torsional  displacement  between  the  mount  and  the  engine  block  will  introduce 
errors  and  hence  are  to  be  minimized.  The  two  previous  generations  of  optical  encoder 
mounts  allowed  the  encoders  to  move  too  much  during  engine  operation.  The  first 
encoder  mount  did  not  provide  adequate  support  for  the  encoder  bearings  [Ref.  7]  The 
second  encoder  mount  was  essentially  an  aluminum  plate  mounted  on  the  end  of  three 
cantilevered  beams.  When  a  proximeter  was  used  to  measure  the  displacement  of  the 
aluminum  plate  with  respect  to  the  engine  block,  it  was  determined  that  the  optical 
encoder  mount  was  moving  approximately  ±0.02  degrees  with  respect  to  the  engine 
block  across  a  wide  variation  of  speeds  and  torques.  This  relative  displacement,  however 
small  it  seems,  corresponds  to  a  20%  error  in  delta  t  readings  with  a  3600  count  optical 
encoder,  and  a  4%  error  in  deltat  readings  from  a  720  count  optical  encoder.  These 
magnitudes  of  error  are  unacceptable  when  changes  in  deltat's  between  different  angular 
positions  can  vary  by  as  little  as  three  microseconds.  A  new  more  rigid  optical  encoder 
mount  was  designed  and  successfully  installed. 

Appendix  B  describes  the  design  of  the  new  optical  encoder  mount  and  optical 
encoder  installation  procedure.  Appendix  F  is  the  procedure  for  establishing  top  dead 
center  (TDC)  relative  to  the  optical  encoder.  All  events  for  data  collection  are  referenced 
to  TDC  of  the  number  one  piston.  Figure  3.2  shows  the  new  optical  encoder  mount. 
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Figure  3.2  New  optical  encoder  mount 


2. 


Pressure  Transducers 


These  pressure  transducers  are  used  to  measure  gas  pressure  inside  the  cylinders 
for  cylinder  pressure  versus  crank  position  data  collection.  Cylinder  pressure  is  a 
necessary  component  in  the  solution  of  the  model  predicted  response 

Reference  air  pressure  is  a  required  input  for  the  Engine  Cycle  Analyzer  software. 
These  pressure  transducers  required  a  reference  pressure  to  properly  measure  the  absolute 
pressure.  This  reference  pressure  is  the  air  box  pressure  or  the  air  pressure  the  cylinder 
experiences  when  the  intake  ports  are  exposed  (see  Figure  3.1).  A  pressure  comparator 
connected  to  the  spark  plug  mounting  hole  in  the  side  of  the  engine  was  used  to  establish  a 
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reference  pressure  versus  engine  speed  curve.  Figure  3.3  is  the  measured  pressure  in  the 
air  box  versus  engine  speed  curve.  Appendix  C  describes  the  procedure  used  to  establish 
the  curve. 


500 


Reference  Air  Box  Pressure  vs.  Engine  Speed 


1000  1500 

Engine  Speed  (RPM) 


2000 


2500 


Figure  3.3  Reference  pressure  versus  speed  curve 


3.  Modulation  Domain  Analyzer 

A  Hewlett  Packard  533 10A  Modulation  Domain  Analyzer  (MDA)  was  used  to 
measure  the  time  intervals  between  successive  windows  of  the  optical  encoder  and  hence 
determine  6{t),  6{t),  and  0{t) .  This  device  measures  frequency  or  time  interval 
(deltat)  between  successive  events  versus  time  (the  modulation  domain).  The  optical 
encoder  installed  on  the  engine  provides  an  output  signal  of  a  repeating  five  volt  TTL 
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square  wave  with  each  rise  corresponding  to  a  known  angular  position  (0. 1  degrees  for  a 
3600  count  optical  encoder  and  0.5  degrees  for  a  720  count  optical  encoder).  With  the 
32,000  bit  extended  memory  option,  eight  revolutions  of  crankshaft  information  using  the 
3,600  count  optical  encoder  and  44  revolutions  of  information  using  the  720  count  optical 
encoder  (8  revs  x  3600  events  per  rev  =  28,800  events  and  44  revs  x  720  events  per  rev  = 
31,680  events)  could  be  recorded  and  stored.  This  information  is  used  to  compare 
measured  angular  position  to  model  predicted  angular  position  in  Section  IV.  C. 

The  procedures  for  measuring  deltat's  from  the  optical  encoder,  accessing  this 
data  from  the  memory  of  the  MDA,  and  storing  it  to  facilitate  passage  to  software  such  as 
MATLAB  is  discussed  in  Appendix  D. 

4.  Data  Acquisition  Procedure 

Figure  3.4  is  a  schematic  of  the  test  engine  and  instrumentation.  Appendix  E 
describes  the  data  acquisition  procedure. 
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Figure  3.4  Test  engine  and  instrumentation 

C.         ENGINE  CYCLE  ANALYZER 

The  Engine  Cycle  Analyzer  (ECA)  manufactured  by  SuperFlow  is  a  PC  based  data 
acquisition  system  which  allows  pressure  versus  angular  position  of  the  crankshaft  data  to 
be  measured  and  stored.  During  data  acquisition,  the  ECA  was  used  to  record  and  store 
pressure  data  from  all  three  cylinders.  The  ECA  has  the  ability  to  monitor  1  to  999 
revolutions  of  data  from  the  engine  and  store  the  ensemble  average  of  this  data  to 
memory.  When  using  the  3,600  count  optical  encoder,  8  revolutions  of  data  were 
recorded  and  when  using  the  720  count  optical  encoder,  44  revolutions  of  data  were 
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recorded.  These  number  of  cycles  were  chosen  to  match  the  maximum  number  of  cycles 
that  the  MDA  could  store.  Revolutions  of  recorded  data  corresponded  to  the  data 
collection  capability  of  the  MDA.  Appendix  E  discusses  operation  of  the  EC  A  and  setup 
files. 
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IV.  RESULTS 

A.         MEASURED  SHAFT  SPEED  VARIATIONS 

To  obtain  a  sufficient  data  base,  time  resolved  cylinder  pressures  and  near 
instantaneous  shaft  speed  were  acquired  for  all  combination  of  speeds  and  torques  shown 
in  Table  4.1.  These  represent  low,  medium,  and  high  torque  and  low,  medium,  and  high 
speed.  For  each  condition  (e.g.  1,000  RPM,  135  FT-LBS)  data  were  acquired  with  both 
the  720  count  and  the  3,600  count  optical  encoder. 

Table  4. 1  Data  runs 


SPEED 
(RPM) 

MEASURED  LOAD  TORQUE 
(FT-LBS) 

1000 

22 

80 

135 

1500 

22 

135 

160 

2000 

50 

135 

200 

Only  the  result  of  the  1,000  rpm  and  135  ft-lb  run  utilizing  the  720  count  optical 
encoder  was  chosen  for  analysis.  Raw  data  from  the  Modulation  Domain  Analyzer  was 
phase  locked  and  ensemble  averaged  to  represent  one  typical  cycle  (0  to  360  degrees  after 
TDC).  Appendix  G  is  the  MATLAB  code  which  phase  locks  and  ensemble  averages  the 
raw  data.  Figure  4.1  is  a  plot  of  the  phase  locked,  ensemble  average  data  from  the  MDA 
for  44  cycles  taken  at  1,000  RPM  and  135  FT-LBS  torque. 
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Figure  4. 1  MDA  phase  locked,  ensemble  average  of  raw  data  taken  at  1,000  rpm 
and  135  ft-lbs. 


Figure  4. 1  shows  on  the  ordinate  the  time  between  successive  windows  (deltat) 
for  the  optical  encoder  versus  time.  This  corresponds  to  360  degrees  of  data.  If  the 
crankshaft  were  turning  with  a  constant  angular  velocity,  the  data  shown  in  Figure  4.1 
would  be  a  horizontal  line.  Any  value  below  the  average  deltat  indicates  a  speed-up  of 
the  crankshaft  (it  takes  a  shorter  time  for  the  next  window  to  arrive).  Any  value  above  the 
average  deltat  indicates  a  slow  down  of  the  crankshaft  (it  takes  more  time  for  the  next 
window  to  arrive).  The  average  angular  velocity  for  the  1,000  rpm,  135  ft-lb  run  is 
107.268  rad/sec.  Appendix  G  is  the  MATLAB  code  that  computes  average  angular 
velocity.  The  average  time  between  pulses  is  81 .35  micro-seconds. 
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These  data  can  be  presented  in  any  of  several  ways.  The  reciprocal  of  deltat 
yields  the  shaft  speed.  However,  these  data  tend  to  be  quite  noisy.  Therefore,  the  data 
will  be  processed  and  presented  using  the  deviation  in  angle  from  a  constant  speed  shaft. 

With  the  ensemble  averaged  deltat's  and  the  average  angular  velocity  computed, 
the  measured  response  at  the  free  end  of  the  crankshaft  can  be  calculated  and  plotted.  The 
equation  to  calculate  is: 

&  measured  =  9 krvwn  ~  ^ cycle  * l  measured    where:    0  blown  =  angular  position  of  optical  encoder 

window  relative  to  TDC 
®  cycle  =  average  angular  velocity  of  crankshaft 

(4.1)  *  measured  =  ^z to  arrival  of  optical  encoder 

window 

Plotting  the  angle  of  the  crankshaft  as  opposed  to  the  angular  velocity  response  reduces 

noise.  Taking  the  integral  of  the  angular  velocity  produces  results  that  are  less  sensitive  to 

noise.    If  the  crankshaft  was  operating  at  a  constant  angular  velocity,  the  response  plot 

would  be  a  horizontal  line.   Deviation  from  constant  velocity  indicates  that  the  crankshaft 

angular  position  is  leading  or  lagging  a  constant  velocity  phasor.     Appendix  G  is  the 

MATLAB  code  that  computes  measured  crankshaft  position  deviation  at  the  free  end  of 

the  crankshaft.    Figure  4.2  is  a  plot  of  measured  crankshaft  position  deviation  from  a 

constant  speed  shaft. 
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Figure  4.2  Measured  Response  of  crankshaft  at  free  end 

B.         TORSIONAL  VIBRATION  MODEL  PREDICTED  RESULTS 

Utilizing  the  measured  pressures  and  calculated  inertia/stiffness  from  the  model, 
the  equations  of  motion  were  solved  numerically  utilizing  a  fourth  order  Runge-Kutta 
method.  The  initial  conditions  were:    0{t  =  0)  =  0  rad  and  0(t  =  0)-  108.9  rad/sec  for  all 

degrees-of-freedom.  The  material  damping  was  chosen  as  0.002  and  viscous  damping 
0.02.  The  non-linear  reciprocating  inertia  terms  were  updated  every  0.5  degrees  or  720 
times  per  revolution.  The  system  of  12  first  order  differential  equations  was  integrated 
from  tinitiai  =  0  to  tfmai  =  0.058  seconds.  There  were  720  linearly  spaced  time  steps  between 
tinitiai  and  tfmai.    These  time  steps  were  chosen  to  match  the  ECA  data  updates.    At  the 
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completion  of  each  time  step,  the  response  of  9i  was  recorded.  Appendix  G  is  the 
MATLAB  code  that  computes  predicted  response.  Figure  4.3  is  a  plot  of  predicted 
crankshaft  position  deviation  from  a  constant  speed  shaft. 
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Figure  4.3  Model  predicted  response  of  crankshaft  at  free  end 

C.         COMPARISON  OF  PREDICTED  TO  MEASURED  RESPONSE 

Figure  4.4  shows  a  comparison  of  predicted  and  measured  response  versus  time. 
Several  similarities  are  obvious.  First,  the  slope  reversals  occur  at  nearly  the  same  phase 
and  time.  Secondly,  the  slopes  are  generally  the  same.  Measured  and  predicted  responses 
show  good  qualitative  agreement. 
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Figure  4.4  Measured  and  predicted  responses  of  crankshaft  at  free  end 

To  better  understand  these  curves  it  is  necessary  to  consider  the  driving  torque. 
Figure  2.5  shows  that  the  maximum  gas  torque  for  the  number  one  cylinder  occurs  20 
degrees  after  TDC.  With  the  crankshaft  rotating  at  107.3  rad/sec  (1024  rpm),  20  degrees 
after  TDC  of  number  one  cylinder  occur  approximately  0.0033  seconds  into  the  cycle. 
The  system  speeds  up  (positive  slope)  until  the  applied  torque  is  zero  at  approximately  80 
degrees  (0.013  seconds)  after  TDC.  When  the  applied  torque  is  zero  the  slope  of  the 
predicted  response  reverses.  This  pattern  continues  for  the  firings  of  cylinders  three  and 
two.  All  three  distinct  increases  in  angular  velocity  of  the  crankshaft  associated  with 
cylinders  firing  are  clearly  present. 
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The  measured  response  with  all  of  the  stated  similarities  still  has  one  large 
difference.  At  the  beginning  of  the  cycle,  the  slope  decreases  and  continues  this  pattern 
until  0.025  seconds  after  TDC.  After  0.025  seconds,  the  angular  velocity  increases  until 
the  cycle  repeats.  Examining  the  measured  response  curve  indicates  that  the  applied 
torque  into  the  system  is  lower  at  the  beginning  of  the  cycle  and  then  increases  toward  the 
end  of  the  cycle  allowing  repeatability.  Figure  2.5,  the  individual  cylinder  gas  torque 
curves,  show  low  gas  torque  for  the  number  one  cylinder,  slightly  increased  gas  torque  for 
number  three  cylinder,  and  finally  the  greatest  gas  torque  contribution  by  number  two 
cylinder.  This  would  appear  to  explain  the  shape  of  the  measured  response  curve,  low 
pressure  in  number  one  cylinder. 

To  validate  the  measurement,  pressure  transducers  were  interchanged  between 
cylinders.  The  same  pressure  curves  were  obtained. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.         SUMMARY 

A  new  rigid  optical  encoder  mount  was  designed  and  installed  on  the  engine.  A 
torsional  vibration  model  was  developed  which  included  the  crankshaft,  major 
reciprocating  components,  and  external  loads.  Equations  of  motion  describing  the 
dynamics  of  the  system  were  derived.  The  Modulation  Domain  Analyzer  and  Engine 
Cycle  Analyzer  were  used  to  gather  time  resolved  information  from  the  engine.  Baseline 
data  were  obtained  at  eighteen  different  speed  and  torque  combinations  utilizing  a  720  and 
3,600  count  optical  encoder.  The  predicted  response  and  the  measured  response  were 
compared. 

C.        CONCLUSIONS 

The  predicted  and  measured  response  at  the  free  end  of  the  crankshaft  show  good 
qualitative  agreement.  The  predicted  response  shows  the  expected  increase  and  decrease 
in  instantaneous  angular  velocity  as  a  result  of  cylinders  firing.  The  torsional  vibration 
model  with  non-linear  inertia  and  linear  stiffness  and  damping  appears  to  describe  the 
dynamics  of  the  diesel  engine  well.  Specifically,  the  six  degrees-of-freedom  is  sufficient  to 
capture  the  measured  torsional  dynamics.  Measured  response  appears  to  be  different  due 
to  actual  low  cylinder  pressure  in  number  one  cylinder. 
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C.         RECOMMENDATIONS 

The  applied  torque  to  the  shaft  should  be  calculated  using  a  parasitic  torque,  where 
a  constant  magnitude  is  subtracted  from  each  cylinder  instead  of  multiplying  each  gas 
torque  by  the  same  mechanical  efficiency.  This  should  greatly  improve  the  agreement 
between  the  model  and  the  measurements. 

Other  contributions  of  inertia  from  the  camshaft,  counter-balance  shafting,  gears, 
and  exhaust  valves  should  be  explored  for  significance. 

More  model  validation  can  be  obtained  by  instrumenting  the  flywheel  with  a 
proximeter.  This  adds  instantaneous  angular  velocity  information  from  the  fifth  degree  of 
freedom.  The  torsional  twisting  between  the  first  and  fifth  degrees  of  freedom  may  also 
hold  information  of  interest. 

The  low  pressure  in  the  number  one  cylinder  should  be  investigated  and  corrected. 
Once  corrected,  another  1,000  rpm,  135  ft-lb  torque  data  run  should  be  made.  The 
measured  response  curve  can  be  compared  to  the  measured  response  from  this  research. 
It  is  expected  that  the  results  will  be  similar  to  the  model  predicted  response. 

The  next  step  should  be  pressure  prediction  based  upon  instantaneous  shaft  speed 
measurement.  The  measured  response  curve  could  also  be  used  to  localize  cylinder  and 
engine  faults. 
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APPENDIX  A 
INERTIA  AND  STIFFNESS  FOR  THE  TORSIONAL  ENGINE  MODEL 

To  calculate  the  equivalent  inertia  and  stiffness  of  the  model,  Reference  1  was 
utilized.  In  general,  the  complex  system  must  be  converted  into  a  simpler  system  of 
equivalent  masses  and  equivalent  shafts.  The  manufacturer's  design  drawings  [Ref.  9]  of 
the  crankshaft  were  utilized  to  obtain  the  geometry  and  material  properties. 

MASS  POLAR  MOMENT  OF  INERTIA  OF  CRANKSHAFT  AND  RECIPROCATING 
COMPONENTS: 

1 .  Each  equivalent  mass  polar  moment  of  inertia  of  the  crankthrows  include  all  inertia 
contributions  from  the  mid-point  of  a  journal  bearing  to  their  left  and  to  the  their  right, 
Figures  2.3  and  2.4  germane.  The  inertia  from  the  mid-point  of  the  front  journal  bearing 
surface  to  the  end  of  the  crankshaft  is  Ji.  It  includes  the  inertia  contribution  of  the  optical 
encoder  crankshaft  adaptor.  The  inertia  from  the  mid-point  of  the  rear  journal  bearing 
surface  to  the  back  of  the  crankshaft  including  the  flywheel  is  J5.  The  equivalent  inertia  of 
the  water  brake  dynamometer  is  J6.  The  value  of  J6  was  obtained  from  the  manufacture  of 
the  water  brake  [Ref.  13].  The  reciprocating  inertia  from  the  pistons  and  connecting  rods 
is  discussed  later. 

2.  In  general,  to  compute  the  equivalent  inertia  of  the  system,  the  crankshaft  must  be 
broken  down  into  simpler  components.  The  following  are  the  component  inertias  which 
must    be    summed    to    obtain   the    equivalent    inertia    at    each    degree    of  freedom. 
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a.  Crankshaft  journals  (solid): 


W  =  -*Df*A*p  (lbf) 


Kl=^-   (inz) 
8 

i      W*K*     nuc-  2s 

J= (lbf-in-sec  ) 

g 
where:  W  =  weight  of  journal  in  lbf 

K  =  radius  of  gyration  of  journal  in  inches 

Di  =  diameter  of  journal  in  inches 

A  =  length  of  journal  in  inches 

p  =  specific  weight  of  material  =  0.283  lbf  in3 

g  =  386  in/sec2 


b.  Crankpins  (solid): 

(lbf-in-sec2) 


D]*B 


1.735 


D2 

^-  +  R2 
8 


where:  D2  =  diameter  of  crankpin  in  inches 
B  =  length  of  crankpin  in  inches 
R  =  distance  of  center  of  gravity  of  crankpin 

from  the  crankshaft  axis  of  rotation,  in  inches 

c.  Crankwebs: 

J=yW*K  (M_m_SQC2) 

g 

where:  W  =  total  weight  of  component  in  lbf 

K  =  radius  of  gyration  in  inches 

g  =  386  in/sec2 
note:     Each  crankweb  must  be  broken  down  into  small  geometric 

shapes.  The  parallel  axis  theorem  is  utilized  to  transfer 

the  equivalent  inertia  to  the  crankshaft  axis  of  rotation. 

Table  10.13  [Ref  1]  contains  radii  of  gyration  (K)  and 

weights  (W)  for  28  geometric  shapes. 

d.  Reciprocating  inertia: 

WR2  n,c.         2. 

Jrec  = (lbf-in-sec  ) 

2*g[l-cos(2<9)] 

where:  W  =  weight  of  piston  and  connecting  rod  =10.5  lbf 

R  =  distance  of  center  of  gravity  of  crankpin 

from  the  crankshaft  axis  of  rotation,  in  inches 

g  =  386  in/sec2 

0  =  angular  position  of  crankthrow  after  TDC 
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EQUIVALENT  TORSIONAL  STIFFNESS  OF  CRANKSHAFT: 

1 .  Each  equivalent  torsional  stiffness  was  calculated  from  one  degree  of  freedom  to 
the  next.  Torsional  stiffness  between  crankthrows  is  from  the  mid-point  of  the  crankpin  of 
one  cylinder  to  the  mid-point  of  the  crankpin  of  the  next  cylinder. 

2.  In  general,  to  compute  the  equivalent  torsional  stiffness  of  the  system,  the 
crankshaft  must  be  broken  down  into  equivalent  lengths  of  shafting.  Component 
equivalent  lengths  must  be  summed  to  get  equivalent  lengths  between  each  degree  of 
freedom. 

a.  For  journals  (solid),  crankpins  (solid),  and  crankwebs: 


Le=  D* 


A  +  0A*DA      \B  +  0.4 * D2  \      {R-02(D,+D2) 


one  journal 


one  pin 


two  crankwebs 


(in) 


where:  D  =  diameter  of  equivalent  shaft  in  inches 
A  =  length  of  journal  in  inches 
Di  =  diameter  of  crank  journal  bearing  in  inches 
D2  =  diameter  of  crankpin 
R  =  distance  of  center  of  gravity  of  crankpin 

from  the  crankshaft  axis  of  rotation,  in  inches 
T  =  thickness  of  crankweb  in  inches 
W  =  width  of  crankweb  in  inches 

(lbf-in/rad) 


K  = 


32*4 
where:  G  =  modulus  of  rigidity  =  12.0  X  106  (lbffin2  ) 


b.  For  solid  shafts: 

tt*D4*G 


K 


32*1. 


(lbf-in/rad) 


c.  To  calculate  the  equivalent  stiffness  from  crankpin  to  crankpin,  calculate  the 
stiffness  as  shown  above  (Kmain)    Main  indicates  main  journal  bearing  surface  to 
main  journal  bearing  surface.  Kp^o-pin  =  (2*Kmam)"1  +  (2*Kmain)"1 
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d.  Figure  11.5  [Ref.  1]  contains  some  typical  crankshaft  elements  and  equivalent 
lengths. 
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APPENDIX  B 

OPTICAL  ENCODER  AND  OPTICAL  ENCODER  MOUNT  INSTALLATION 

PROCEDURE 


In  order  to  obtain  accurate  readings  from  the  optical  encoder,  the  following 
procedure  should  be  followed  when  installing  the  mount  and  or  changing  optical  encoders: 
Installing  mount: 

1 .  Place  the  crankshaft  to  optical  encoder  stainless  steel  adapter  on  the  end  of  the 
crankshaft  and  firmly  seat.  Rotate  the  adapter  until  the  two  encoder  shaft  set  screws  are 
pointing  straight  up  and  down  or  vertical  to  the  deck.  Tighten  the  three  set  screws 
securely  against  the  crankshaft  with  an  alien  wrench. 

2.  Place  the  smaller  of  the  two  aluminum  optical  encoder  mounts  against  the  engine 
block  and  then  place  the  larger  one  on  top  oriented  such  that  the  encoder  will  mate. 
Obtain  the  optical  encoder  mount  bolts  and  hand  tighten  the  bolts  until  the  mounts  are 
flush  with  the  block.  Figure  3.2  shows  this  orientation. 

3.  Obtain  the  optical  encoder  alignment  jig.  Place  it  into  the  countersunk  encoder 
receptor.  The  purpose  of  the  jig  is  to  align  the  mount  so  the  crankshaft  rotates 
concentrically  with  the  optical  encoder  shaft.  Carefully  tighten  down  all  three  optical 
encoder  mounting  bolts.  Ensure  that  the  jig  is  able  to  move  in  and  out  of  the  encoder 
adapter  hole  and  the  mounting  lip  without  binding.  It  may  be  necessary  to  shim  the 
encoder  mount  If  all  three  bolts  are  securely  fastened  and  the  jig  moves  in  and  out  freely 
from  the  receptor,  the  optical  encoder  mount  is  aligned. 
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Installing  optical  encoder: 

1 .  Place  optical  encoder  into  the  receptor  and  install  the  three  small  bolts  which  retain 
it.  Tighten  the  bolts  with  an  alien  wrench. 

2.  On  the  bottom  of  the  optical  encoder  there  is  a  window  which  allows  access  to  the 
encoder  shaft  set  screws.  Tighten  the  one  set  screw  securely  against  the  optical  encoder 
shaft  with  an  alien  wrench.  It  may  be  necessary  to  utilize  a  flashlight  and  mirror  to 
accomplish  this. 

3.  Remove  the  access  plate  over  the  flywheel.  Rotate  the  engine  crankshaft 
clockwise  with  a  pry  bar  180  degrees.  This  allows  the  other  the  set  screw  to  be  accessed. 

4.  Tighten  the  remaining  set  screw  securely  against  the  optical  encoder  shaft  with  an 
alien  wrench. 

The  optical  encoder  mount  and  optical  encoder  are  now  aligned.  Ensure  that  TDC 
relative  to  the  optical  encoder  is  established  prior  to  engine  operation,  Appendix  G 
germane. 
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APPENDIX  C 
ESTABLISHING  REFERENCE  PRESSURE 

In  order  to  obtain  the  most  accurate  cylinder  pressure  readings,  the  ECA  requires  a 
reference  pressure  input.  This  appendix  describes  the  procedure  for  establishing  reference 
pressure  versus  speed.  It  is  obtained  at  normal  operating  temperature  of  the  engine. 

An  MKS  Instruments  Inc.  Type  698  differential  high-accuracy  pressure  transducer 
in  conjunction  with  the  MKS  Instruments  Inc.  Signal  Contitioner  Type  270B  was  used  to 
measure  engine  air  box  pressure  relative  to  atmospheric  pressure.  The  range  of  the 
differential  pressure  transducer  was  1,000  torr  (20  psi),  sufficient  to  measure  the  air  box 
pressure.  The  Air  box  pressure  was  obtained  by  removing  a  cold  weather  spark  plug  on 
the  side  of  the  engine  and  using  the  threaded  hole  for  a  air  box  pressure  tap.  The  other 
side  of  the  differential  pressure  transducer  was  open  to  the  atmosphere.  The  engine  was 
operated  from  700  to  2,500  rpm  in  increments  of  200  rpm.  At  every  increment  of  speed, 
the  air  box  pressure  was  recorded.  A  plot  of  reference  air  box  pressure  versus  engine 
speed  was  constructed,  Figure  3.3  germane.  When  measuring  cylinder  firing  pressures, 
the  reference  air  box  pressure  for  the  speed  of  engine  operation  should  be  input  into  the 
ECA  set  up  file. 

The  following  picture  shows  the  set  up  of  the  pressure  transducer,  signal 
conditioner,  and  associated  piping  to  the  engine. 
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Figure  C.  1  Instrumentation  utilized  for  establishing  reference  air  box  pressure 
versus  engine  speed 
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APPENDIX  D 
MODULATION  DOMAIN  ANALYZER  OPERATING  PROCEDURE 

This  appendix  describes  how  to  configure  and  operate  the  modulation  domain 
analyzer  (MDA)  to  record  time  resolved  information  from  the  optical  encoder  during 
engine  operation.  The  procedure  is: 

1 .  Ensure  both  the  MDA  and  data  acquisition  computer  are  turned  off. 

2.  Connect  one  end  of  the  IEEE-488  cable  (National  Instruments  763061-02  Rev  C, 
Type-X2,  Length  2. 1  meter)  to  the  HPIB  on  the  back  side  of  the  MDA.  Connect 
the  other  end  of  this  cable  to  the  National  Instruments  GPIB  board  interface  on 
back  of  the  data  acquisition  computer. 

3.  Connect  the  TDC  output  on  front  of  the  EC  A  sensor  interface  to  the  EXT  Arm 
input  on  front  of  the  MDA  with  a  BNC  connector. 

4.  Connect  the  CA  output  on  front  of  the  EC  A  sensor  interface  to  the  Channel  A 
input  on  front  of  the  MDA  with  a  BNC  connector. 

5.  Turn  the  power  to  the  MDA  on. 

6.  Depress  the  'function  and  input'  menu  key. 

7.  Select  the  'time  interval'  soft  key. 

8.  Select  'common'  A  and  B  line  soft  key. 

9.  Select  'rising  edge  A  and  Rising  edge  B'  soft  key. 

10.  Depress  the  'trigger'  menu  key. 

1 1 .  Select  'triggered'  soft  key. 

12.  Select  'mean'  and  '  1/mean'  on  the  entry  panel  to  display  the  mean  deltat  and 
mean  frequency. 

13.  Start  the  engine  in  accordance  with  Appendix  A  [Ref.  1]. 

14.  Depress  the  'auto  scale'  menu  key. 

15.  The  MDA  is  now  displaying  deltat  on  the  y-axis  and  time  on  the  x-axis.  Continue 
with  step  16  to  make  the  MDA  ready  for  remote  operation  and  data  acquisition. 

16.  Turn  on  the  power  to  the  data  acquisition  computer  and  go  to  the  HP  Applications 
icon  and  enter. 

17.  For  remote  operation,  depress  the  'utility'  menu  key  on  the  MDA. 

1 7.  Select  'addressed'  under  the  HP-IB/Print  soft  key. 

18.  Select  address  number  5  (five). 

19.  The  MDA  is  now  in  remote  operation  mode.  In  order  to  print  the  screen  or 
operate  the  MDA  locally,  the  'TALK  ONLY'  mode  must  be  enabled. 
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The  following  section  discusses  how  to  access  the  memory  of  the  MDA,  record 
the  measured  deltat's,  and  store  this  information  with  the  data  acquisition  computer. 


1 .  With  the  MDA  properly  configured,  select  the  HP  Applications  icon  and  go  to  the 
HP53305A  Phase  Analyzer  software.  Connect  'CA'  lead  to  'B'  input  on  MDA. 

2.  Under  the  setup  menu,  select  instrument  setup  and  configure  the  following: 

a.  Input  configure:  B,  +  slope,  DC  coupling,  and  50  ohm  impedance. 

b.  Arming:  EXT,  TTL,  and  +  slope. 

c.  Interface:  National  GPE3,  0  board  index,  and  address  5. 

4.  Under  the  setup  menu,  select  measure/view  setup  and  configure  the  following: 

a.  Parameter:  Frequency  and  Phase. 

b.  Number  of  samples:  Enter  the  number  of  sample  points  desired  for  the  run  but 
do  not  exceed  32,000.  Remember,  the  720  count  optical  encoder  is  720  data 
samples/revolution  and  the  3,600  optical  encoder  is  3,600  samples/revolution. 

c.  Measurement  time:  Enter  the  approximate  expected  deltat's  expected  for  the 

speed  you  are  taking  data  at.  (i.e.  at  lOOOrpm,  with  the  720  optical  encoder, 

>     Q?  ■>     ■               a            min      *         rev         *60sec 
enter  83.3  microseconds  or * * . 

l,000rev     HOwindows      min 

d.  Number  of  averages:  one. 

c.  Pacing:  auto. 

d.  Input:  B. 

e.  Clock/carrier  frequency:  calculated. 

f.  No  filter. 

5.  Under  the  setup  menu,  select  phase  units  then  select  'unit  interval'. 

6.  The  MDA  is  now  ready  to  collect  TDC  triggered  data  from  the  optical  encoder. 

7.  When  ready  to  start  collecting  data,  go  to  the  measure  menu  and  select  'phase 
deviation'.  As  soon  as  you  release  the  mouse  button,  the  MDA  records  data 

8.  Once  the  data  is  recorded,  the  MDA  will  ask  you  what  to  name  the  plot.  This 
is  unimportant,  enter  anything. 

9.  Save  the  data  to  the  hard  drive  or  zip  drive  as  a  M-File. 

10.  With  MATLAB,  remove  the  header  and  footer  information  which  leaves  two 
columns  of  numbers.  The  length  of  the  array  is  the  number  of  sample  points  you 
previously  selected  during  setup  of  the  software. 

1 1 .  The  data  file  is  now  ready  for  use  by  MATLAB  code,  Appendix  H  germane. 
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Figure  D.  1  MDA  and  high  speed  data  acquisition  computer.  An  oscilloscope  sets 
on  top  of  the  MDA. 
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APPENDIX  E 
DATA  ACQUISITION  PROCEDURE 

This  appendix  describes  the  data  acquisition  procedure  for  comparing  model 
predicted  response  to  measured  response  at  the  free  end  of  the  crankshaft  utilizing  the  720 
count  optical  encoder.  The  procedure  is: 


1 .  Align  the  engine  for  operation  and  turn  the  charge  amplifiers  on.  Ensure  the 
proper  scale  and  sensitivity  is  set  on  charge  amplifiers.  This  information  is  on  the 
Kistler  pressure  transducer  calibration  sheets. 

2.  Setup  the  ECA  setup  file  used  for  data  acquisition  in  accordance  with  Appendix  F. 

3.  Align  the  MDA  and  data  acquisition  computer  in  accordance  with  Appendix  D. 

4.  Increase  the  engine  and  load  to  the  desired  set  point. 

5.  On  the  ECA  select  data  taking  and  statistics  under  the  main  menu. 

6.  Select  720DD3 53 . SET  under  the  data  taking  and  statistics  menu.  720DD3 53 . SET 
is  the  preset  for  monitoring  and  recording  all  three  cylinder  pressures  utilizing  the 
720  count  optical  encoder.  DD353PRE.SET  is  the  setup  file  for  monitoring  and 
recording  all  three  cylinder  pressures  utilizing  the  3,600  count  optical  encoder. 

7.  Enter  the  correct  reference  pressure  from  Figure  3.3  into  the  J.  1  block  and  enter 
the  number  of  cycles  to  be  samples  (should  correspond  with  the  MDA). 

8.  Place  the  charge  amplifiers  in  the  'operate'  mode. 

9.  Enter  'G'  for  graphic  monitoring  and  ensure  three  pressure  versus  theta  curves 
appear. 

10.  The  ECA  is  now  ready  to  record  data. 

1 1 .  The  ECA  will  temporarily  turn  on  and  off  the  output  crank  angle  and  TDC  signals 
while  it  is  acquiring  data.  Therefore,  the  MDA  data  must  be  taken  first. 

12.  Record  the  speed  and  torque  of  the  engine  by  depressing  'Record  Data'  on  the 
SuperFlow  Dynamometer  Control  Console. 

1 3 .  Take  the  MDA  measured  data  in  accordance  with  Appendix  D. 

14.  Immediately  after  the  MDA  data  acquisition  data  is  taken,  depress  "D"  for  data 
taking  on  the  ECA  keyboard.  Follow  on  screen  instructions. 

15.  Store  each  individual  cylinder  pressure  to  the  hard  drive. 

16.  Data  acquisition  is  complete.  Steps  12  through  15  above  require  practice  and 
should  be  accomplished  as  fast  as  possible  to  avoid  any  deviation  in  speed  and 
torque  from  the  engine. 
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Figure  E.  1  EC  A,  EC  A  sensor  interface,  and  SF-901  Dynamometer  computer. 


Figure  E.2  SuperFlow  Dynamometer  Control  Console 
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APPENDIX  F 

ESTABLISHING  TOP  DEAD  CENTER  RELATIVE  TO  THE  OPTICAL 

ENCODER 


In  order  to  get  accurate  readings  from  the  optical  encoder,  the  following 
procedure  should  be  followed  to  correctly  position  Top  Dead  Center  (TDC)  of  the 
number  one  cylinder  relative  to  the  optical  encoder.  TDC  is  the  reference  against  which 
all  angular  measurements  are  made.  It  is  the  point  in  the  cycle  of  a  cylinder  when  the 
volume  is  the  least.  The  number  one  cylinder  is  the  cylinder  for  this  reference.  The 
procedure  outlined  in  Appendix  G  of  Reference  1  is  unchanged  with  the  following 
exceptions. 

1 .  To  rotate  the  engine,  remove  the  flywheel  access  cover  on  the  back  of  the  engine. 

2.  A  pry  bar,  as  the  one  depicted,  can  be  utilized  to  rotate  the  engine. 

3.  Once  mechanical  TDC  of  the  number  one  cylinder  is  established,  loosen  the  three 
retaining  bolts  that  secure  the  optical  encoder  to  the  mount. 

4.  Rotate  the  optical  encoder  until  the  TDC  signal  on  the  EC  A  sensor  interface 
illuminates. 

5.  Carefully  tighten  the  retaining  bolts  on  the  optical  encoder  base. 
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Figure  F.  1  Flywheel  access  cover  removed  with  pry  bar  in  place. 
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APPENDIX  G 
MATLAB  CODE  USED  FOR  CYCLE  ANALYSIS 

This  Appendix  contains  all  the  different  MATLAB  M-files  used  for  this  research. 
Data  from  both  the  ECA  and  MDA  require  that  the  footer  and  header  information  be 
removed  prior  to  loading  into  an  M-file.  The  ECA  records  data  pressure  from  -180 
degrees  to  180  degrees  relative  to  TDC.  Before  loading  the  ECA  data  into  an  M-file  for 
analysis,  it  must  first  be  rearranged  and  saved  from  0  degrees  to  360  degrees  after  TDC. 
This  procedure  must  be  followed  for  each  cylinder  pressure  data. 

The  following  code  can  be  utilized  to: 

a.  Plot  pressure  versus  crank  angle  for  all  three  cylinders. 

b.  Plot  gas  torque  versus  crank  angle  for  all  three  cylinders. 

c.  Compute  average  gas  torque  contribution  from  each  cylinder 

d.  Compute  the  mechanical  efficiency  based  upon  a  torque  balance. 

e.  Compute  and  plot  the  phase  lock,  ensemble  average  of  MDA  raw  data. 

f.  Compute  the  average  angular  velocity  through  each  revolution  from  MDA  data. 

g.  Solve  the  equations  of  motion  of  the  system  and  plot  a  predicted  response  at 
the  free  end  of  the  crankshaft. 

h.  Compute  and  plot  the  measured  response  at  the  free  end  of  the  crankshaft 
utilizing  MDA  phase  lock,  ensemble  average  data. 

I.  Define  the  equations  of  motion  as  a  system  of  12  first  order  differential 
equations. 
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%  This  code  computes  and  plots  the  Individual  cylinder  pressures 

%  Referenced  to  Nr.  1  cylinder  TDC 

%  Filename:  cylpress.m 

load  dl5krl35.m    %  ECA  cylinder  #1  pressure  data 

load  el5krl35.m    %  ECA  cylinder  #2  pressure  data 

load  fl5krl35.m    %  ECA  cylinder  #3  pressure  data 

% 

%PLOTS  PRESSURE  OF  EACH  INDIVIDUAL  CYLINDER  VS  DEGREES  AFTER  TDC  OF  NR1  CYL 

% 

theta=l:  1:720; 

plot(theta/2.(dl5krl35(theta)*14J),,x,,theta/2,(el5krl35(theta)*14.5).,o,.theta/2,(fl5krl35(theta)*14.5),,- 

-') 

title('Individual  Cylinder  Pressure  referenced  to  TDC  of  Nr.  1  cylinder') 

xlabel('Degrees  after  TDC  of  Nr.  1  Cylinder') 

ylabel('Pressure  (PSI)') 

legend('Nr  1  cyl  pressure','Nr  2  cyl  pressure', 'Nr  3  cyl  pressure') 

axis([-inf3610  1250]) 


%  This  code  computes  the  individual  torque  contribution  of  each 

%  individual  cylinder  referenced  to  TDC  of  Nr.  1  Cylinder  and  the 

%  recipricating  inertia.  It  can  also  be  manipulated  to  plot 

%  pressure  vs.  theta  curves 

%  Filename:  torque. m 

% 

load  dlkrl35.m   %  ECA  cylinder  #1  pressure  data 

load  elkrl35.m    %  ECA  cylinder  #2  pressure  data 

load  flkrl35.m   %  ECA  cylinder  #3  pressure  data 

r=2.25;W=10.5;g=386;  %constants 

for  theta=  1:720 
Tlcyl(theta)=(dlkrl35(theta)*14.5*r*sin(theta*.008726646)*(pi*rA2))/12;  %FT-LB 

T2cyl(theta)=(elkrl35(theta)*14.5*r*sin((theta*.008726646)-(4*pi/3))*(pi*rA2))/12;%FT-LB 
T3cyl(theta)=(flkrl35(theta)*14.5*r*sin((theta*.008726646)-(2*pi/3))*(pi*rA2))/12;%FT-LB 
j2rec(theta)=(W*rA2/(2*g))*(l-cos(2*theta*008726646));  %lb*in*secA2 

j3rec(theta)=(W*rA2/(2*g))*(l-cos(2*(theta*.008726646)-(2*pi/3)));    %lb*in*secA2 
j4rec(theta)=(W*rA2/(2*g))*(l-cos(2*(theta*.008726646)-(4*pi/3)));    %lb*in*secA2 

%   plot(theW2,(Tlcyl+T2Cyl+T3cylVkX') 

%   titleCNet  Torque  Input  Referenced  to  TDC  of  Nr.  1  Cylinder') 

%   xlabel('Degrees  after  TDC  of  Nr.  1  Cylinder') 

%   ylabel('Gas  Torque  (FT-LB)') 

%   hold  on 

end 

figure(l) 

theta=linspace(  1 ,720,720); 

plot(theta/2.Tlcyl,'X',theta/2,T2cyl,,0',theta/2,T3cyl,'.') 

title('Individual  Cylinder  Torque  input  Referenced  to  TDC  of  Nr.  1  Cylinder') 

xlabelCGas  Torque  (FT-LB)') 

ylabel('Degrees  after  TDC  of  NR.  1  Cylinder') 
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legend('Nr  1  Cyl  TorqueVNr  2  Cyl  TorqueVNr  3  Cyl  Torque') 

figure(2) 

plot(theta/2j2rec.,X'.theW2.j3rec.'0'.theta/2.j4rec/.') 

title('Individual  Reciprocating  Crankthrow  Inertia  Referenced  to  TDC  of  Nr.  1  Cylinder') 

xlabel('Inertia  (lb*in*secA2') 

ylabel('Degrees  after  TDC  of  Nr.  1  Cylinder') 

legend('Nr  1  Cyl  Reciprocating  Inertia','Nr  2  Cyl  Reciprocating  Inertia'.'Nr  3  Cyl  Reciprocating  Inertia') 


% 

%  THIS  SECTION  COMPUTES  THE  AVERAGE  TORQUE  CONTRIBUTION  OF  EACH  CYLINDER 

% 

tic 

x=linspace(0,(2*pi),720); 

T  lcyl_avg=trapz(x,T  lcyl)/(2*pi) 

T2cyl_avg=trapz(x,T2cyl)/(2*pi) 

T3  cyl_avg=trapz(x,T3cyl)/(2  *pi ) 

AvgTorque  input=(T  1  cyl_avg+T2cyl_avg+T3  cyl_avg)/3 

Torque_out_to_Torque_in=134/(Tlcyl_avg+T2cyl_avg+T3cyl_avg) 


%  Jim  Hudson    Thesis  file 

%  720  OE,  SAMPLE  RUN  1000  RPM.  135  FT-LB  TORQUE 

%  16  NOV  1997  Filename:  ensemble,  m 

%  3 1680  sample  points  from  MDA  or  44  revs  of  INFO 

%  THIS  CODE  PHASE  LOCKS  AND  ENSEMBLE  AVERAGES 

%  THE  OUTPUT  REPRESENTS  ONE  REVOLUTION  OF  INFO 

load  7 1  kr  1 3  5 .  m  %  This  is  where  the  3 1 680  times  are  located 

t=X7 lkrl35(:,l);  %  Assigns  vector  "t"  to  each  time  sample 

for  b=l:3 1680 

tt(b)=t(b+ 1  )-t(b);      %  Computes  the  delT's  and  assigns 
end  %  them  to  Vector  "tt" 

tt=(reshape(tt, 720,44))';  %  Phase  locks  and  adds 
tttt=sum(tt)/44;  %  Averages  the  phase  lock 

position=linspace(0.(2*pi),720); 
timer=0;  %  Initialize  a  counter 

for  e=  1:720 

time=[e;t(e+l);tttt(e)]; 

plot(timer.tttt(e).'kX')        %PLOTS  TIME  VS.  DEL_T 

timer=timer+tttt(e); 

hold  on 
end 

gtext('1000  RPM  &  135  Ft-Lbs,  720  Optical  Encoder*) 
title('MDA  DATA  PHASE  LOCK  ENSEMBLE  AVERAGED') 
ylabel('Delta_t  btwn  windows  (Seconds)') 
xlabel(Time  of  sample,  after  TDC  (SECONDS)') 


%  This  code  computes  the  average  angular  velocity  through  each 
%  of  the  44  revolutions  of  of  MDA  data.  Additionally,  it  gives 
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%  the  average  angular  velocity  the  MDA  saw. 

%  filename:  omegacyc.m 

% 

load  71krl35.m  %  This  is  where  the  MDA  raw  data  is  stored 

for  a=l:3 1681 

t(a)=X71krl35(a); 
end 
for  n=  1:44 

x(n)=720*n; 
end 

omega_zero=(2*pi)/(t(720)-t(  1)); 
plot(0,omega_zero.'kX') 
hold  on 

sum=omega_zero; 
for  c=  1:43 

omega_avg(c)=(2  *pi)/(t(x(c+ 1  ))-t(x(c)+ 1 )) ; 

sum=sum+omega_avg(c); 
end 

c=linspace(  1,43,43); 

average=sum/44         %  This  is  the  average  angular  velocity 
plot(c,omega_avg,'kX') 

title(' Average  Angular  Velocity  Through  Each  of  the  44  Revolutions') 
xlabel('Cycle  Number") 
ylabelC Average  Angular  Velocity  (rad/sec)') 


%  THIS  CODE  PLOTS  CONCURRENTLY  THE  MEASURED  AND  PREDICTED  RESPONSES 

%  for  a  720  count  optical  encoder.  The  measured  and  predicted 

%  response  can  be  plotted  independently  by  simply  commenting  out 

%  one  or  the  other.      Filename:  measpred.m 

%  3 1680  sample  points  from  MDA  are  located  in  file:  71KR135.M 

%  ECA  cylinder  pressure  information  is  in  files: 

%  dlkrl35.m  cylinder  #1  pressure  data 

%  elkrl35.m  cylinder  #2  pressure  data 

%  flkrl35.m  cylinder  #3  pressure  data 

% 

%  This  first  section  plots  the  measured  response 

% 

load  71krl35.m  %  THIS  IS  WHERE  THE  3 1680  TIMES  ARE  LOCATED 

t=X71krl35(:,  1);         %  ASSIGNS  VECTOR  "t"  TO  EACH  TIME  SAMPLE 

for  b=  1:3 1680 

tt(b)=t(b+ 1  )-t(b);     %  COMPUTES  THE  DEL_T'S  AND  ASSIGNS  THEM  TO 
end  %  VECTOR  "tt" 

tt=(reshape(tt,720,44))';%  This  phase  locks  one  cycle 
tttt=sum(tt)/44;         %  This  ensemble  averages  the  phases 
position=linspace(0,(2*pi),720);  %  The  known  positions  of  the  OE.  windows 
timer=0; 
for  e=  1:720 

plot(timer,(position(e)-(107.278*timer)),'rO')  %107.278  is  omega_bar 

timer=timer+tttt(e) ; 

hold  on 
end 
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gtext('Measured  Response') 
hold  on 
% 

%  This  section  is  a  SIX  SECOND  ORDER  SIMULTANEOUS  EQUATION  ODE  SOLVER 
%  This  section  calls  "eqns.m"  which  defines  the  system  of  2nd  order  ode's 
% 

global  Tload  plcyl  p2cyl  p3cyl  Tlcyl  T2cyl  T3cyl ... 
J2recj3recj4recjl  J2j3j4j5j6kl  k2  k3  k4... 

k5  cl2  c23  c34  c45  c56  c2  c3  c4  %share  the  variables  with  other  m-files 
ic=[0  0  0  0  0  0  108.909  108.909... 

108.909  108.909  108.909  108.909];  %define  IC's  [x(t=0)  &  xdot(t=o)] 
load  dlkrl35.ni         %  cylinder  #1  pressure  data  from  the  ECA 
load  elkrl35.m  %  cylinder  #2  pressure  data  from  the  ECA 

load  flkrl35.m         %  cylinder  #3  pressure  data  from  the  ECA 
Tload=  1608;  %  lb-in  ( 1 34  FT-LB  as  measured  by  the  dynamometer) 

z=linspace(0,.057692066.721);  %  divide  .058  sees  into  720  evenly 

%  spaced  time  increments 
r=2.25;     %  Crankshaft  eccentricity  (in) 
W=  10.5;    %  The  weight  of  the  reciprocating  components 

%  (piston  and  connecting  rod)  (lbf) 
g=386;      %  The  acceleration  of  gravity  (in/secA2) 
mech_eff=.6175;  %  This  is  the  scale  factor  for  the  Gas  Torque  in 

%mech_eff=Tload/(Tlavg+T2avg+T3avg) 

%  M-file  Torque. m  computes  the  gas  torques 
omega_bar=109.9;  %This  is  the  mean  speed  of  the  crankshaft  used 

%by  the  plot  command  for  "Response  Deviation" 
j  1=0291;     %lb*in*secA2/rad 
j2=.2730;     %lb*in*secA2/rad 
j3=.1710;     %lb*in*secA2/rad 
j4=.2730;     %lb*in*secA2/rad 
j5=6.194;     %lb*in*secA2/rad 
j6=.2800;     %lb*in*secA2/rad 
kl=3.470e6;  %lb*in/rad 
k2=9.551e6;  %lb*in/rad 
k3=9.551e6;  %lb*in/rad 
k4=13.495e6;  %lb*in/rad 
k5=1.304e6;  %lb*in/rad 
cl2=002;     %lb*in*sec/rad 
c23=.002;     %lb*in*sec/rad 
c34=.002;     %lb*in*sec/rad 
c45=002;     %lb*in*sec/rad 
c56=.002;     %lb*in*sec/rad 
c2=.02;       %lb*in*sec/rad 
c3=02;       %lb*in*sec/rad 
c4=.02;       %lb*in*sec/rad 

x=ic;         %Initialize  the  "x"  vector  to  compute  Tcyl's  and  Jrec's 
for  theta=  1:720 

Tlcyl=mech_eff*dlkrl35(theta)*14.5*r*sin(theta*.008726646)*(pi*rA2);%lb-in 
T2cyl=mech_eff*elkrl35(theta)*14.5*r*sin((theta*.008726646)-(4*pi/3))*(pi*rA2);%lb-in 
T3cyl=mech_eff*flkrl35(theta)*14.5*r*sin((theta*.008726646)-(2*pi/3))*(pi*rA2);%lb-in 
j2rec=(W*rA2/(2*g))*(l-cos(2*theta*008726646));  %lb*in*secA2 

j3rec=(W*rA2/(2*g))*(l-cos(2*(theta*.008726646)-(2*pi/3)));    %lb*in*secA2 
j4rec=(W*rA2/(2*g))*(l-cos(2*(theta*.008726646)-(4*pi/3)));    %lb*in*secA2 
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[T.x]=ode45('eqns\z(theta).z(theta+ 1  ).ic); 

ic=x(length(T),:);  %  REINITALIZE  THE  IC'S  FROM  THE  PREVIOUS  ITERATION 
plot(T(length(T)).x(length(T),l)-(omega_bar*T(length(T))).,bX,) 
end 

title('Crank  Position  (PREDICTED/MEASURED)  Deviation  from  Mean") 
xlabel('Time  (sec)') 
ylabel('Angular  Position  (rad)') 

axis([0  .0573  -.027  .01])  %  Change  this  for  the  appropriate  desired  view 
gtextCPredicted  Response') 


%  This  code  defines  the  12  first  order  equations  to  be  differentiated. 

%  FILENAME:  EQNS.M 

function  xdot=firsteqn(t.x,Tload) 

global  Tload  Tlcyl  T2cyl  T3cyl  j2rec  j3rec  j4rec  jl  j2  j3... 

j4  j5  j6  kl  k2  k3  k4  k5  cl2  c23  c34  c45  c56  c2  c3  c4 
% 

%  These  are  12  first  order  equations  which  define  the 
%  original  6  second  order  equations  of  motion! ! ! ! 
% 

xdo(l)=x(7); 
xdo(2)=x(8); 
xdo(3)=x(9); 
xdo(4)=x(10); 
xdo(5)=x(ll); 
xdo(6)=x(12); 

xdo(7)=((cl2*x(8))-(cl2*x(7))+(kl*x(2))-(kl*x(l)))/jl; 
xdo(8)=(Tlcyl+(c23*x(9))-((c23+cl2+c2)*x(8))+(cl2*x(7))+(k2*x(3))... 

-((kl+k2)*x(2))+(kl  *x(  1  )))/(j2+j2rec); 
xdo(9)=(T2cyl+(c34*x(  10))-((c34+c23+c3)*x(9))+(c23  *x(8))+(k3  *x(4)). . . 

-((k2+k3)*x(3))+(k2*x(2)))/03+j3rec); 
xdo(10)=(T3cyl+(c45*x(ll))-((c45+c34+c4)*x(10))+(c34*x(9))+(k4*x(5))... 

-((k3+k4)*x(4))+(k3*x(3)))/04+j4rec); 
Xdo(ll)=((c56*x(12))+((c56+c45)*x(ll))+(c45*x(10))+(k5*x(6))... 

-((k4+k5)*x(5))+(k4*x(4)))/j5; 
xdo(  1 2)=(-Tload-(c56*x(  1 2))+(c56*x(  1 1  ))-(k5*x(6))+(k5*x(5)))/j6; 
xdot=xdo';  %  THIS  IS  THE  VECTOR  DEFINING  THE  EQUATIONS  OF  MOTION 
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